Quantum mechanics insights into melatonin and analogs binding to melatonin MT1 and MT2 receptors

Melatonin receptors MT1 and MT2 are G protein-coupled receptors that mediate the effects of melatonin, a hormone involved in circadian rhythms and other physiological functions. Understanding the molecular interactions between these receptors and their ligands is crucial for developing novel therapeutic agents. In this study, we used molecular docking, molecular dynamics simulations, and quantum mechanics calculation to investigate the binding modes and affinities of three ligands: melatonin (MLT), ramelteon (RMT), and 2-phenylmelatonin (2-PMT) with both receptors. Based on the results, we identified key amino acids that contributed to the receptor-ligand interactions, such as Gln181/194, Phe179/192, and Asn162/175, which are conserved in both receptors. Additionally, we described new meaningful interactions with Gly108/Gly121, Val111/Val124, and Val191/Val204. Our results provide insights into receptor-ligand recognition’s structural and energetic determinants and suggest potential strategies for designing more optimized molecules. This study enhances our understanding of receptor-ligand interactions and offers implications for future drug development.

www.nature.com/scientificreports/In this context, in silico studies have been described as an important initial step in the development of new, more potent drugs based on existing ones 54,55 .These studies involve the analysis of protein-ligand interactions using energetic calculations methodologies through classical mechanics (MM), quantum mechanics (QM), or hybrid methods (QM/MM) [56][57][58][59] .This allows for the evaluation of the strengths and weaknesses of each ligand, proposing structural modifications that enhance its affinity for the receptor.
Molecular dynamics (MD) simulations using calculations based on classical mechanics are widely used to obtain optimized structural models of receptor-ligand complexes that are not available in the Protein Data Bank (PDB), especially those obtained through in silico molecular docking 60,61 .From the ensemble of conformations generated in MD, protein-ligand interaction calculations can be obtained with relatively low computational cost using methodologies such as Quantum-Mechanical/Molecular-Mechanical Generalized Born Surface Area (QM/ MM-GBSA), which is based on the use of semi-empirical functional for quantum mechanics calculations, and classical field force for molecular mechanics calculations.
In this case, semi-empirical methods are reduced-order Hamiltonians in which some elements are replaced by empirical parameters adjusted to experimental or ab initio data.Examples include the PM3, PM6, and AM1 models, among others 62,63 .Always aiming for increased accuracy and reduced computational cost, semi-empirical functionals are applied to only a part of the complex (usually defined by the ligand distance), while the remaining is calculated using classical force fields.Despite being a methodology used in drug development studies 64 , it is known that its accuracy is still lower, especially when compared to purely ab initio methods such as Density Functional Theory (DFT) 65 .
The DFT emerged as a quantum functional to address the limitations of semi-empirical approaches and other quantum methodologies, such as Hartree-Fock (HF), enabling the study of biological molecules 66,67 .Despite the advantages of DFT, it remains computationally demanding, requiring the use of high-performance clusters and other cost-reduction techniques.One such approach is molecular fragmentation with conjugated caps (MFCC), which involves splitting the system into smaller molecules 71 .
In this aspect, the present work investigates the interaction of two melatonin receptors (MT 1 and MT 2 ) complexes with RMT, 2-PMT, and MLT.For this purpose, the crystallographic data of the structural complexes and complexes obtained from docking simulations were used as input geometrical data in the calculation.The quantum binding energy calculations were performed using the DFT formalism DFT in association with the MFCC method.This approach allows us to determine the main contributions of amino acids to their affinity with melatonin receptors.

Molecular docking and molecular dynamics simulation
To obtain the MT 1 -MLT and MT 2 -MLT complexes, molecular docking was performed using AutoDock Vina 68 .The distribution of Vina scores (in kcal/mol) from 1000 rounds of docking can be observed in Fig. 1A, and  www.nature.com/scientificreports/summary of the obtained values is presented in Table 1.The Vina scores for the MT 1 -MLT complex ranged from − 7.607 kcal/mol to − 8.235 kcal/mol, with a mean (SD) value of − 7.832 kcal/mol (0.08753).Similarly, for the MT 2 -MLT complex, the Vina scores ranged from − 7.420 kcal/mol to − 8.021 kcal/mol, with a mean (SD) of − 7.632 kcal/mol (0.07979).These results indicate that melatonin has a higher affinity for the MT 1 protein compared to MT 2 .However, it is important to note that the primary objective of the exhaustive Vina docking is to rank the complexes based on Vina scores.Furthermore, as depicted in Fig. 1B, the best MLT conformations bind to the same pocket as observed in experimental structures 23 , further confirming the reliability of the docking results.Therefore, the obtained data strongly suggest that the docking results are suitable for proceeding to MD simulations.The main purpose of the MD simulations is to adjust the conformation of the protein and ligand in a solvent environment.The Root Mean Square Deviation (RMSD) is a metric that evaluates the difference between 3D structures based on atomic distances.A higher RMSD value indicates a greater difference in structures.In this study, we compared the structures along the trajectory with the initial frame.The RMSD of the MT 1 -MLT complex can be seen in Fig. 2A.It is noteworthy that all three replicates exhibited similar behavior throughout the 200 ns trajectory, with an RMSD value of approximately 0.5 nm.Additionally, the Root Mean Square Fluctuation (RMSF) analysis (Fig. 2C), which evaluates the average residue fluctuation, showed no significant differences among the three replicates.It is expected to observe higher fluctuation in the N-and C-termini of the protein.However, in this case, only higher residue fluctuation was observed in the C-terminal region of the protein.
Unlike the MT 1 -MLT complex, the three replicates of MT 2 -MLT exhibited distinct behaviors along the MD trajectory.Among them, the RMSD (Fig. 2B) of replicate 2 was more stable.Like MT 1 -MLT, the RMSD value remained around 0.5 nm.In contrast, replicate 3 showed a higher but stable RMSD (0.65 nm) until 125 ns into the simulation.Afterward, the RMSD increased to values higher than 1.1 nm, and after 150 ns, stability was observed near this RMSD value.The heatmap plot of RMSD per residue over time (Supplementary Figure S1C) revealed an increasing fluctuation of residues attached in the MT 2 N-terminal, which corresponds to the apocytochrome BRIL (UniProt P0ABE7) from Escherichia coli, after 125 ns.Additionally, around 60 ns, an increased RMSD in residues of MT 2 C-terminal corresponding to another fusion protein, rubredoxin (Rub, UniProt P00268).Notably, these residues are distant from the ligand binding pocket, suggesting that they have insignificant influence on the ligand binding mode.Finally, replicate 1 showed intermediate fluctuations with values near 0.8 nm until 175 ns, followed by a small decrease to values near 0.70 nm.The fusion protein also observed these variations, as shown in Supplementary Figure S1A.The structure of MT 2 and the fusion proteins can be seen in Supplementary Figure S2.The RMSF analysis of the MT 2 replicas (Fig. 2D) indicates that the fluctuation patterns were quite similar among them, despite the amplitude of replica 3 being larger, which corroborates with the observations made in the RMSD graph.Conversely, the fluctuations of replicas 1 and 2 were alike in terms of both amplitude and fluctuation pattern.The RMSD analysis of the ligand is an interesting means to evaluate if the ligand remained in the binding pocket or if displacement occurred.A high RMSD indicates a modification in the ligand's position.Supplementary Figure S3 shows that MLT in complex with MT 1 and MT 2 remained in the pocket as the RMSD did not show high values.However, MLT in MT 1 showed an RMSD near 0.5 nm, while MLT in MT 2 showed an RMSD near 0.8 nm, suggesting that MLT-MT 2 underwent a higher conformational change during the MD simulation compared to MLT-MT 1 .These results were found to be suitable for performing the QM/MM-GBSA analysis.
The results of QM/MM-GBSA can be observed in Table 2 and Fig. 3.In Table 2, it is possible to observe that the mean of all replicas showed lower energy in the interaction with MLT for MT 1 compared to MT 2 , except for replica 3.This is also true when we compare the lowest value for each replicate.Notably, the QM/MM-GBSA analysis of MT 1 -MLT (Fig. 3A) showed more variation among replicates compared to the MT 2 -MLT complex.However, both complexes exhibited values mostly lower than − 10 kcal/mol.As observed in the Vina docking results, the complex with higher affinity was MT 1 -MLT.The selected frames for MFCC/DFT analysis were frame 1648 (replicate 2, MT 1 ), which presented − 27.57kcal/mol of binding energy, and frame 1871 (replicate 2, MT 2 ) which presented − 24.78 kcal/mol of energy.Now, there are six complexes for which QM calculations will be performed: MT 1 -MLT, MT 2 -MLT, MT 1 -RMT, MT 2 -RMT, MT 1 -2-PMT, and MT 2 -2-PMT.The last four complexes were obtained from the Protein Data Bank (accession codes: 6ME2, 6ME9, 6ME3 and 6ME6 respectively).
The structures of MT 1 -MLT and MT 2 -MLT selected for QM (DFT) calculations were compared with their experimentally resolved counterparts, both in their active and inactive forms.The structure of MT 1 -MLT was compared with structures PBD ID: 7DB6 73 (active, Cryo-EM) and PDB ID: 6ME2 25 (inactive, X-ray) and the complex MT 2 -MLT was overlaid with structures PDB ID: 7VH0 69 (active, Cryo-EM) and 6ME6 23 (inactive, X-ray).As observed in Fig. 4A, a displacement of TM6 in the active state (light blue ribbon) is noted when compared to MT 1 with MLT (light green ribbon).The same is observed in MT 2 -MLT (Fig. 4B). Figure 4C and D present the overlay of MT 1 -MLT and MT 2 -MLT with inactive MT 1 and MT 2 , respectively.It is noted that the displacement of TM6 when comparing the complexes with each other is much smaller, indicating that the complex with MLT is anchored in the inactive form of the melatonin receptor.This is expected since the docking was performed with the protein in this conformation.Therefore, it is concluded that the MD simulation maintained the complex in its inactive state.

MFCC and quantum mechanical calculations
To aid in examining and analyzing binder interactions, we schematically divided the binders into regions, as illustrated in Fig. 5.The carbon atoms are labeled from 1 to 19, and the rings are denoted by letters (A, B, and C).According to the Marvin Sketch analysis, neutral states were predominant (100%) at both pH values examined (7.0 and 7.4).In Fig. 5B, the MT 1 and MT 2 receptors are represented with the RMT and 2-PMT molecules, respectively, highlighting the binding site of the ligands.
Energy calculations and evaluations were performed for each of the six complexes, and convergence criteria were determined based on the variation of total energy (in kcal/mol) and radius r (in Å).The cumulative binding  energy was obtained by summing up all the calculated energies between the residues within the radius and the binder.The analyses were continued until the contributions of amino acids from each radius were no longer significantly different from the total interaction energy found after each successive radius, with a threshold of less than 10%.Figure 6 illustrates the change in total interaction energy with increasing radius.The convergence criteria are met between r = 7 Å and 8 Å for all complexes and dielectric constants (ε = 10 and ε = 40).However, calculations were conducted within r = 10 Å to ensure the comprehensive evaluation of crucial residues.As a result, for the MT 1 -RMT and MT 2 -RMT complexes, 96 and 105 residues were found to interact within r = 10 Å, respectively.In terms of affinity, the total interaction energy when ε = 40 (ε = 10) for MT 1 -RMT was − 43.08 kcal/mol (− 43.84 kcal/mol) and for MT 2 -RMT − 48.75 kcal/mol (− 50.07 kcal/mol), suggesting that the RMT drug has a slightly higher affinity for MT 2 -RMT.In the case of MT 1 -2-PMT and MT 2 -2-PMT, binding affinity was assessed for 104 and 108 residues, respectively.The total energies for these residues were − 60.Based on the data above, it is observed that the increasing order of affinity of the MT 1 (MT 2 ) receptor is: RMT < MLT < 2-PMT (MLT < RMT < 2-PMT), which suggests that 2-PMT acts better on both receptors compared to RMT and MLT.Moreover, all energy values show what is expected for ε = 10 compared to ε = 40, where the first shows lower values due to increased medium permittivity.Furthermore, it is observed that the pattern for ε = 10 and ε = 40 remained consistent in the plots, which strongly suggests that the DFT calculations were carried out correctly.Therefore, from now on, we will present the study with energy results focused on the dielectric constant ε = 40.

MT 1 interaction analysis
The individual energetic contribution of each amino acid was evaluated to understand the most significant residues in the interaction with MT 1 , as well as the ligand regions where these interactions occurred.In Fig. 7, the calculated energy of the key amino acids that contributed (both positively and negatively) to the interaction energy between the complexes can be observed.For interactions within the MT 1 -RMT complex (Fig. 7, green bar), the ten most important residues in decreasing order of affinity (in kcal/mol) are the following: Phe179 (− 7. Seven common amino acid residues were observed between them: Phe179, Met107, Gly108, Val111, Val191, Val159, and Ile112.These residues are then suggested as key interacting residues, and main interactions with them will be evaluated in detail. Phe179 was the most important residue for interaction in all three MT 1 complexes.It formed a pi-alkyl interaction with atom [ii(C8)] of the RMT molecule (Fig. 8A) and two pi-pi interactions with rings [A,B] of 2-PMT (Fig. 8B) and MLT (Fig. 8C).Interestingly, Phe179 showed high affinity for ligands (between − 7.04 kcal to − 9.21 kcal/mol) regardless of whether they interacted with its aromatic ring or aliphatic carbon.Val159 interacted with atoms from the aromatic region of RMT and 2-PMT (Fig. 8A-B) through dipole-dipole interactions.In contrast, it formed a non-conventional hydrogen bond with MLT (Fig. 8C), which resulted in a higher affinity interaction (− 4.30 kcal/mol) than the other molecules.Ile112 exhibited dipole-dipole interactions with all molecules.However, the interaction energies were lower for RMT (Fig. 8A) and 2-PMT (Fig. 8B) (− 1.93 kcal/ mol and − 2.01 kcal/mol, respectively) than for MLT (Fig. 8C) (− 3.31 kcal/mol).This could be attributed to the distance factor, as Ile112 interacted with the same shared region of MLT (2.8 Å) and 2-PMT (3.6 Å).For RMT, the interaction involved a hydrogen atom bonded to a carbon atom of the aromatic ring [ii(C14)H], which may have affected the charge distribution and attractive forces.Met107 was another important key residue for interaction with all ligands.It showed lower energy values for RMT (− 3.83 kcal/mol) (Fig. 8D) and 2-PMT (− 3.70 kcal/mol) (Fig. 8E) than for MLT (− 1.93 kcal/mol) (Fig. 8F).It formed a pi-sulfur interaction with ring [B] of RMT (Fig. 8D), a dipole-dipole interaction and a non-conventional hydrogen bond with atom [ii(N)H] of 2-PMT (Fig. 8E), and a dipole-dipole interaction with atom [ii(N)H] of MLT (Fig. 8F).RMT did not meet the non-conventional hydrogen bond parameters.Val111 performed pi-alkyl interactions with aromatic rings of RMT (Fig. 8D) and 2-PMT (Fig. 8E), which resulted in higher affinity than its dipole-dipole interaction with MLT (Fig. 8F).Gly108 exhibited considerable affinity for all ligands.It formed dipole-dipole interactions with atoms [ii(C15)H] of RMT (Fig. 8D), [ii(C10) and ii(C11)] of 2-PMT (Fig. 8E), and [ii(N)] of MLT (Fig. 8F).The latter was also a hydrogen bond.Finally, Val191 showed  MT 2 interaction analysis Among the residues of MT 2 -RMT complex (Fig. 9, green bar), the ten residues that most contributed to the binding affinity in decreasing order (in kcal/mol) were: Phe192 (− 9.As observed for MT 1 , Phe192 of MT 2 structure had the highest affinity for RMT (− 9.95 kcal/mol) and 2-PMT (− 9.22 kcal/mol).It formed two pi-pi interactions with two aromatic rings of each ligand (Fig. 10A-B).In the MT 2 -MLT complex, Phe192 had a lower energy value (− 1.79 kcal/mol) and a dipole-dipole interaction with atom [i(C1)H] (Fig. 10C).Val124 performed a pi-alkyl interaction with RMT (− 2.56 kcal/mol), a pi-sigma interaction  with the aromatic ring of 2-PMT (− 4.12 kcal/mol), and a dipole-dipole interaction with atom [ii(C5)] of MLT (− 2.84 kcal/mol) (Fig. 10A-C).Gly121 showed affinity energies lower than − 3.0 kcal/mol for RMT and 2-PMT (− 3.22 kcal/mol and − 3.19 kcal/mol, respectively).It formed a dipole-dipole interaction with atom [ii(C9)H] of RMT and a non-conventional hydrogen bond with atom [ii(N)] of 2-PMT (Fig. 10A,B).For MLT (− 1.09 kcal/ mol), it formed a non-conventional hydrogen bond with atom [i(O)] (Fig. 10C).

Alanine scanning
The alanine scanning study is an exciting approach to assessing the significance of amino acids in protein-ligand complexes.The PremPLI server, a straightforward tool based on machine learning, predicts the mutation's effect solely using 3D structural information from the complex.In this case, the evaluated interactions involved the amino acids Gly108/Gly121, Val111/Val124, and Val191/Val204 across all ligands (RMT, 2-PMT, and MLT).The analysis of amino acids Gln181/194, Phe179/192, and Asn162/175 was restricted to those ranked with low energy according to QM/DFT calculations.
As seen in Table 3, the mutation with the most significant impact on interaction reduction was F179A/F192A.The highest quantum calculations energy for Phe192 was associated with the MLT ligand, and according to PremPLI, it also exhibited the most negligible impact on the mutation.Notably, among the amino acids Gly108/ Gly121, Val111/Val124, and Val191/Val204, the G108A/G121A mutation had the most substantial impact on interaction, according to PremPLI.However, it is observed that in all complexes, the mutations resulted in positive ∆∆G values, indicating reduced affinity.Despite PremPLI's results aligning with QM calculations, further robust analyses are essential to confirm the importance of amino acids Gly108/Gly121, Val111/Val124, and Val191/Val204 in ligand recognition and interaction with MT 1 and MT 2 receptors.

Discussion
The search for new therapeutic strategies for the treatment of insomnia has become extremely important, mainly due to the various adverse effects caused by the use of BZDs and n-BZDs drugs in the general population.Several agonists and analogs of MLT have been studied for years, such as RMT, 2-iodomelatonin, agomelatonin, and 2-PMT 23 .However, all these drugs were developed before the three-dimensional resolution of the MT 1 and MT 2 receptors, using the ligand-based drug discovery (LBDD) strategy.
Thanks to the recent advances in the resolution of the structures of the MT 1 and MT 2 proteins, it is now possible to meticulously evaluate the interactions between the receptors and their ligands.Therefore, this study was developed to fill two gaps related to the melatonin receptors: to obtain a model with MLT and its receptors and to evaluate the biochemical interactions with RMT, 2-PMT, and MLT to propose a structure-based drug discovery (SBDD) strategy, which is reported to be more complex and potentially better than LBDD 70 .The www.nature.com/scientificreports/docking analysis allowed the orientation of MLT in the MT 1 and MT 2 receptors to be predicted in silico, since the structures of the experimental complexes are not yet available.Thus, it was possible to determine how MLT interacts with its receptors.The MD simulation provided important information regarding the stability of MLT in solution when bound to the MT 1 and MT 2 receptors.This is crucial because a molecule that does not stabilize in the active site of the protein with the conformational changes observed in solution does not become a suitable therapeutic target for use.
In this study, we can observe that MLT has good stability when complexed with its targets.In the docking and MD simulation study, it was possible to obtain complexes of MT 1 -MLT and MT 2 -MLT in their inactive state, as well as the other structures analyzed in this study.This comparison was feasible because the TM6 domain in the active form exhibits a displacement relative to the inactive form, which is well described by previous studies with experimentally resolved structures 69 .
In the QM/DFT analysis using the MFCC method, it was observed that the increasing order of affinity for MT 1 was RMT < MLT < 2-PMT, which suggests that 2-PMT is the best ligand for MT 1 activation and RMT is the weaker.In the case of MT 2 this order was: MLT < RMT < 2-PMT, which suggests again 2-PMT as the best binder, however, in this case, MLT was ranked as the weaker MT 2 ligand.All ligands showed higher affinity for MT 1 compared to MT 2 , except the RMT binder.The 2-PMT and MLT results corroborate with in vitro experimental data where the affinity of 2-PMT and MLT for the receptors was evaluated.In these studies, it is possible to observe that the affinity of 2-PMT is higher than the affinity of MLT 71 for both receptors, and both molecules have a higher affinity for MT 1 than for MT 2 [71][72][73] .The RMT results do not agree with the experimental data.Further analysis should be performed to evaluate if this is a limitation of the calculation method or the experimental crystallographic data.However, most of the results shown here are in agreement with in vitro tests.
In the literature, the importance of three amino acids for protein interaction and activation is well described.These residues are Gln181/194, Phe179/192, and Asn162/175 23 .Undoubtedly, the most crucial key residue observed in this study was Phe179/192, which was found to be important in the interaction of all ligands with both proteins.The Gln181/Gln194 was found to be important for RMT (both proteins) and 2-PMT (MT 1 ) interaction, and Asn162/Asn175 was found to be important for RMT (both proteins) and 2-PMT (MT 2 ).Hence, the importance of Gln181/Gln194 and Asn162/175 depends on binder and protein, which is not true for Phe179/192.In addition to these amino acids, three other conserve amino acids-Gly108/Gly121, Val111/Val124, and Val191/ Val204-were shown to play important roles in the interaction with MT 1 and MT 2 for all binders.Beyond Gln181/194, Phe179/192, and Asn162/175 amino acids, only Thr178 and Met107 were recently observed in other studies involving binding with this protein 73 .Thus, we have potential amino acids crucial for protein-ligand interactions with both receptors.This could lead to more optimized molecule development based on the SBBD approach since key residues for the interaction of these molecules were identified, as well as the regions of the melatonin analogs that favored an increased affinity for the receptors.This is an important, although initial, step towards solving the issue of the short half-life of melatonin receptor ligands.
In conclusion, recent advances in deciphering the structures of MT 1 and MT 2 proteins have allowed for a detailed exploration of interactions with their binders.Molecular docking and dynamics simulations were employed to generate MT 1 -MLT and MT 2 -MLT complexes, revealing stronger binding to MT 1 than MT 2 .Notably, 2-PMT displayed higher affinity for both receptors compared to MLT, aligning with experimental findings.The analysis highlighted the importance of specific amino acids, such as Gln181/194, Phe179/192, and Asn162/175, in protein-ligand interactions, providing insights into potential therapeutic strategies.Also, it described new important interactions with Gly108/Gly121, Val111/Val124, and Val191/Val204.This study enhances our understanding of receptor-ligand interactions and offers implications for future drug development.All the analyses of the interactions described here and the key regions of the receptors and ligands involved in the interactions are extremely relevant information for developing more effective drugs for treating sleep disorders.However, it is known that analyses involving drug optimization and mutation analysis are necessary to confirm the research data and fill these gaps.Future studies should consider the virtual screening of molecules with the important regions of the ligands (especially the aromatic rings) to find new ligands with greater affinity and optimize existing molecules by adding and/or removing groups that favor/disfavor interaction affinity.

Molecular docking and molecular dynamics simulation
Since there is no crystallized structure of the MT 1 and MT 2 proteins with MLT, the first step of the study involved obtaining the complexed structures through molecular docking using the AutoDock Vina program (referred to as Vina here) 68 .Initially, the three-dimensional (3D) structure of MLT was obtained from the PubChem website (www.pubch em.ncbi.nlm.nih.gov/compound/896) in SDF format.The molecule's protonation at pH = 7.0 and 7.4 (following the protein experimental pH) was verified using the MarvinSketch code version 17.24 (Marvin Beans Suite -ChemAxon, www.chema xon.com).Subsequently, the molecule was converted to PDBQT format using the OpenBabel 74 server.
The MT 1 structure was obtained from the Protein Data Bank (www.rcsb.org).Four MT 1 crystals were compared for resolution, with codes 6ME2, 6ME3, 6ME4, and 6ME5.The crystal with the lowest resolution (6ME2 25 -2.8 Å), which suggests higher quality, was selected for MLT docking.The same procedure was followed for the MT 2 structure, where structures 6ME6, 6ME7, 6ME8, and 6ME9 were compared.The 6ME6 23 structure was chosen due to its lower resolution (2.8 Å) among MT 2 structures.
Before converting to the PDBQT format, the MT 1 and MT 2 proteins underwent a cleaning step, chain adjustments, and minimization on the Discovery Studios server.In this case, artifacts from the crystallization process were removed, as well as the ligands RMT and 2-PMT, bound to MT 1 and MT 2 , respectively.Additionally, ing side chains and hydrogen atoms were added.The protonations of the MT 1 and MT 2 proteins were evaluated on the online PropKa server for pH 7.0 and 7.4, respectively, in accordance with the crystallization experiment.After these modifications, the protein backbone's conformations were restrained, and an energy minimization (EM) step was performed using the CHARMm (Chemistry at Harvard Molecular Mechanics) force field.The EM utilized convergence tolerances of 10 −5 kcal/mol for total energy change, 10 −3 kcal/mol for the mean square root of the RMS gradient, and 10 −5 Å for the maximum atomic displacement, employing the Smart Minimizer algorithm.For the conversion of large molecules from PDB to PDBQT format, AutoDock Tools 75 was used.
The same binding sites of RMT and 2-PMT were considered for melatonin anchoring.For each system, the melatonin molecule was docked 1000 times, each time generating 10 binding conformation modes.The best conformation (lower Vina score) of each system was selected to proceed to MD simulations.
MD simulations were conducted to optimize the binder conformation within the binding pocket.Independent triplicate simulations were performed using GROMACS 2022 software 81 .The ligand parameters for the MD simulations were obtained using the ACPYPE server (www.bio2b yte.be/ acpyp e/) 76 with Gasteiger as the charge method and GAFF2 as the force field.The server provided all the necessary topology and parameter files for MD simulations using GROMACS software.The force field selected for protein was the Amberff99SB-ILDN.Six MD simulations were carried out for the protein-ligand complex, three for each MT 1 and MT 2 systems.
For each system, a cubic box was utilized with the TIP3P water model extended 12 Å away from solute atoms, and Cl -ions were added to neutralize it.Two rounds of energy minimization were executed to adjust unfavorable contacts in the initial structure.The first minimization step involves a maximum of 20000 steps or until the maximum force on any atom is reduced to below 50 kJ/mol/nm.The steepest descent algorithm was employed with protein restraint to focus on solvent relaxation.The second minimization step, without protein restraint, was performed in flexible water using the same steepest descent algorithm, and the maximum steps were increased to 10,000 or until the force on any atom fell below 250 kJ/mol/nm.The system's pressure and temperature were adjusted to 1 atm and 310 K, in two separate 100 ps steps, referred to as the NVT ensemble (temperature setting) and NPT ensemble (pressure setting).The modified Berendsen 77 and Parrinello-Rahman 78 algorithms were applied to control the system temperature and pressure, respectively.Throughout both steps, hydrogen bonds were constrained using the LINCS algorithm 79 , and positional restraints were applied to the protein to stabilize the solvent around the solute.
Long-range electrostatic interactions were computed using the Particle Mesh Ewald (PME) summation method, employing a non-bonded interaction cut-off of 1 nm.The equations of motion were integrated using the leap-frog algorithm 80 with a time step of 0.2 fs.Before the MD simulation, a small NPT ensemble of 1 ns was conducted without any restrictions on protein position, followed by a production run for 200 ns without restrictions on protein conformation.The MD run produced a total of 2000 protein frames.
The gmx_MMPBSA program 81 was employed to analyze the last 500 frames (50 ns) of the complexes for each MD simulation using hybrid Quantum Mechanics/Molecular Mechanics-Generalized-Born surface area (QM/MM-GBSA) free energy calculation.The QM region was limited to residues near 5 Å from the binder, and the applied semi-empirical functional was PM6-DH + .Solvent molecules and Cl -ions were excluded from the analysis.This method allowed the calculation of interaction energy for various conformations of the biological complex at a relatively low computational cost, enabling the selection of complexes with higher affinities for a more robust and accurate analysis using QM/DFT.Hence, the complexes with lower energy for each system were chosen for further QM (DFT) calculations using the MFCC approach.The MD trajectories were visualized using UCSF Chimera software 82 .Root mean square deviation (RMSD) and fluctuation (RMSF) were calculated using the "gmx" commands of the GROMACS package.All plots were generated using the R language in RStudio 4.1.1(http:// www.rstud io.com/), and protein image representations were created using PyMol 83 .

MFCC and quantum mechanical calculations
The QM calculation is an accurate methodology for studying ligand-protein interactions.However, it has been observed that the computational cost is significantly high for large systems.As the present study aims to analyze protein-ligand interactions, and therefore, is considered a complex biological system for QM calculations, the MFCC methodology developed by Zhang and Zhang 71 was applied to each of these systems to overcome this limitation.The MFCC scheme involves splitting the protein into individual amino acids by breaking the peptide bonds and calculating the interactions between each residue and the ligand separately.The summation of individual amino acid energies provides an approximate binding energy.
The MFCC approach employs "caps" to complete the valence of the amino acids after breaking the peptide bonds at both the N-and C-termini.These caps are composed of amino acid residues that precede and succeed the main amino acid.Additionally, the caps serve the purpose of closely reproducing the amino acid environment.Consequently, it becomes possible to calculate an approximate binding affinity for large systems, such as a protein-ligand complex.
Equation (1) presents the MFCC scheme for calculating the interaction energy (IE(BID/R i )) between the binder (BID) and the amino acid R i , with i denoting the i th amino acid of the protein chain.: The caps, represented as C i−1 and C i+1 , correspond to the neighboring residues covalently bonded to the amine and carboxyl groups of R i , respectively.The first term (E(BID +C i−1 R i C i+1 )) calculates the interaction energy between the binder BID and the main residue (R i ) bonded to the caps (C i-1 and C i+1 ).
However, MFCC aims to assess the individual amino acid contributions to the binding affinity.Thus, the second and third terms are introduced to isolate the contribution of the main residue and eliminate the interaction energy between the binder and the caps.The second term (E(C i−1 R i C i+1 )) represents the energy of the residue (R i ) bonded to the caps (C i-1 and C i+1 ), while the third term (E(BID +C i−1 C i+1 )) indicates the interaction energy between the BID and the caps (C i-1 and C i+1 ).Subtracting the energies calculated in the second and third terms from the first term removes their influence on the energy interaction between BID and R i .
Since the caps' energies were subtracted twice (in the second and third terms), it becomes necessary to add them back in the fourth term (E(C i−1 C i+1 )) to account for their effect on the energy interaction between BID and R i .This step accurately evaluates the individual amino acid's contribution to the overall binding affinity.
For the MFCC involving the MLT molecule, the frames with the lowest energy obtained from the QM/MM-GBSA analysis of MT 1 and MT 2 were selected.No modification was needed since the structures were adjusted previously.As for the 2-PMT molecule, the structures were obtained from the PDB, where the structures of MT 1 and MT 2 were submitted under accession codes 6ME3 25 and 6ME6 23 , respectively.The crystallographic structures of the MT 1 -RMT and MT 2 -RMT complexes were taken from the PDB (PDB codes: 6ME2 25 and 6ME9 23 , respectively).As mentioned earlier, structural adjustments are necessary to address potential limitations of the X-ray technique.Thus, the same cleaning procedure, protonation verification at pH 7.0 (MT 1 ) and 7.4 (MT 2 ) (for both ligands and receptors), the addition of missing side chains, and energy minimization, as performed previously, were applied to these structures.The only difference is that the 2-PMT and MLT molecules were retained in the structures.
After performing fragmentation for each amino acid, the interaction energy between the receptor and the binder was computed using the Gaussian 16 package 84 , which employs the density functional theory (DFT) formalism 85,86 .The simulations were carried out using the generalized gradient approximation (GGA) with the B97D functional, which has been proven to be an efficient and accurate QM method for large systems, particularly when dispersion forces play a significant role 87 .To represent the Kohn-Sham orbitals for all electrons, a small basis set 6-311+G(d,p) with triple split valence (valence triple-zeta), an additional diffuse function (+), and polarization functions (d,p) were applied.The conductor-like polarizable continuum model (CPCM) was utilized to account for solvent effects in the QM calculations, using dielectric constants (ε) of 10 and 40.The constant of 40 is reported as being similar to the crystalline environment 57,58,88 .Meanwhile, the constant of 10 is used here as a control, where a lower constant is expected to result in a higher medium permissivity and lower energy values.Thus, we can ensure that calculations for both constants were performed correctly.Each term of Equation (1) was obtained from the DFT simulations.
It has been established that as the amino acids are increasingly distant from the binder, their contribution to the binding energy with the small molecule is reduced.To eliminate the calculation of unimportant interactions, a convergence analysis of the total binding energy was conducted to limit the number of amino acid residues considered.The analysis incorporated the closest residues from the binding site and excluded the most distant ones.The interaction energy of amino acid residues within imaginary spheres with a pocket radius of r centered at the binder was evaluated, where r = R/2 (for R = 1, 2, 3, 4,..., N n ; with N n being the next natural number of sequences) 88 .The convergence criteria were met when the total energy, while increasing the radius r, did not change by more than 10% compared to the previous radius value. (1)

aFigure 1 .
Figure 1.Melatonin docking in MT 1 and MT 2 .(A) Box plot of Vina score distribution in kcal/mol.The pink box is related to MT 1 -MLT complex and the blue is related to MT 2 -MLT docking.(B) Best docking pose of MLT (yellow stick) bound in MT 1 (pink cartoon) and MT 2 (blue cartoon).

Figure 2 .
Figure 2. RMSD and RMSF analysis of 200 ns three MD trajectory.(A) and (B) RMSD plots in three replicates of MT 1 and MT 2 MD simulations, respectively.(C) and (D) RMSF plots in three replicates of MT 1 and MT 2 by residue number, respectively.

Figure 4 .
Figure 4. Comparison of the MT 1 and MT 2 structures in their active and inactive forms.(A) Comparison of the active MT 1 structures (blue color ribbon) with the MT 1 -MLT complex (green color ribbon).(B) Comparison of the active MT 2 (blue color ribbon) with the MT 2 -MLT complex (green color ribbon).(C) Comparison of the inactive MT 1 structures (pink color ribbon) with the MT 1 -MLT complex (green color ribbon).(D) Comparison of the inactive MT 2 structures (purple color ribbon) with the MT 2 -MLT complex (green color ribbon).

Figure 6 .
Figure 6.Plot representation of total interaction energy of MT 1 and MT 2 complexes for both dieletric constants (ε = 10 and ε = 40) as a function of the binding pocket radius calculated using the DFT/MFCC approach.

Figure 7 .
Figure 7. Graphic panels showing the most important residues for MT 1 interaction with RMT (green bar), 2-PMT (orange bar) and MLT (yellow bar).Also, the region (i, ii or iii) and atom (based on Fig. 5 schematic representation) that interact with each residue at the binding site.

Figure 8 .
Figure 8. Tri-dimensional binding mode of RMT (A and D), 2-PMT (B and E) and MLT (C and F) in MT 1 structure.The ligands are represented in ball and stick (RMT-green, 2-PMT-orange and MLT-yellow and the main MT 1 amino acids in magenta sticks.The interaction types are colored as follows-green: hydrogen bond; cyan: pi-pi interaction; dark blue: non-conventional hydrogen bond; brown: dipole-dipole; red: pi-alkyl; orange: pi-sulfur.

Figure 9 .
Figure 9. Graphic panels showing the most important residues for MT 2 interaction with RMT (green bar), 2-PMT (orange bar) and MLT (yellow bar).Also, the region (i, ii or iii) and atom (based on Fig. 5 schematic representation) that interact with each residue at the binding site.

Table 1 .
Summary of Vina docking results based on 1000 ligand conformations.

Table 3 .
Results of alanine scanning performed by the PremPLI server for the key amino acids indicated by QM/DFT calculations.